Validation of a two-fluid turbulence model in comsol multiphysics for the problem of flow around aerodynamic profiles

The article presents a study of a two-fluid turbulence model in the Comsol Multiphysics software package for the problem of a subsonic flow around the DSMA661 and NACA 4412 airfoils with angles of attack of 0 and 13.87 degrees, respectively. In this paper, the finite element method is used for the numerical implementation of the turbulence equations. To stabilize the discretized equations, stabilization by the Galerkin least squares method was used. The results obtained are compared with the results of other RANS, LES, DES models and experimental data. It is shown that in the case of continuous flow around the DSMA661 airfoil, the results of the two-fluid model are very close to the SST results and are in good agreement with the experimental data. When flowing around the NACA 4412 airfoil, flow separation occurs and a recirculation zone appears. It is shown that in such cases the two-fluid model gives more accurate results than other turbulence models. Implementation of the Comsol Multiphysics software package showed good convergence, stability, and high accuracy of the two-fluid turbulence model.


Abbreviations
Institute of Mechanics and Seismic Resistance of Structures of the Academy of Sciences of the Republic of Uzbekistan, Tashkent, Uzbekistan. 2 Fergana Polytechnic Institute, Fergana, Uzbekistan. 3 Central Aero-Hydrodynamic Institute Named After Professor N.E.Zhukovsky, Zhukovsky, Moscow Region, Russia. 4Moscow Institute of Physics and Technology (MIPT), Dolgoprudny, Moscow Region, Russia.* email: madaliev.me2019@mail.ru1.Using the Custom PDE functions to implement a two-fluid turbulence model in the COMSOL Multiphysics software package.2. Validation of the two-fluid turbulence model and verification of the computational algorithm on a number of simple test problems, such as flows around a flat plate with a zero pressure gradient, a DSMA661 airfoil with an angle of attack of 0 degrees and a NACA 4412 airfoil with an angle of attack of 13.87 degrees.

Two-fluid turbulence model.
The description of this model is presented in several publications by one of the authors of this article [20][21][22] .The main equations for studying the tasks posed are the hydrodynamic equations of the two-fluid model 24 for an incompressible medium In the given system of equations V i is the component of the average flow velocity, ϑ i is the component of the relative velocity, νji is the molar viscosity tensor, p is the pressure, ρ is the density of the medium, ν is the molecular viscosity, K f is the friction coefficient, C s is the coefficient at the Saffman force, def( V ) is the strain rate, determined as: The coefficient of friction is found from the following relation: In this expression, d is the nearest distance to the solid wall, λ max is the real part of the largest root of the characteristic equation det(A − E) = 0, where A is the matrix Constant patterns are C s = 0.2, C 1 = 0.7825, C 2 = 0.306.Consider a two-dimensional stationary solution to system (1).For the finite element method, the application of the standard Galerkin method will lead to a weak form: Here, the equations of system (5) are a weak form of the equation of motion for the averaged and relative velocities, as well as the equation of continuity.Here, v , ṽ , and q are weight functions for the average velocity V , relative velocity ϑ , and pressure p , respectively.
Using the Finite Element Method (FEM) as a discretization method for finding solutions to turbulence models can be quite a challenge.The standard Galerkin problem statement deals with potential sources of numerical instabilities.This is, for example, the case when convection or reaction conditions prevail in the flow 27,28 .It is difficult to find such a stabilization that makes the solution of equations as reliable as possible.A stabilized FEM formulation can be created by adding grid-dependent, consistent, and numerically stabilizing terms to the standard Galerkin method.Various stabilization methods were proposed, many of which are based on the Petrov-Galerkin (PG) approach, which uses a modification of the standard Galerkin weight term.Examples of (1) (3) (4) www.nature.com/scientificreports/popular SG methods are the Streamline-Upwind/Petrov-Galerkin scheme [29][30][31][32][33] (SUPG) and the pressure stabilization scheme/Petrov-Galerkin 34 (SUPG).Over the years, the Petrov-Galerkin method has been developed and gave rise to new, more advanced methods.One of them is the Galerkin Least Squares (GLS) stabilization method.GLS is a general stabilization method applicable to a wide range of problems [35][36][37][38][39] .Its theoretical basis is that the test function should be chosen so as to minimize the squared residual of the equations.By working with matrices and not just with scalar equations, GLS can be formulated for systems of transport equations.By adding the least squares term, we get the Galerkin/least squares formula for the two-fluid equations as shown below: where R GLS and RGLS are given as: Parameters τ and τ for stationary problems are determined by the following expression: And for pressure, parameter τ p is determined by the following expression: where Re h is the Reynolds number of the element:

Solution method and boundary conditions
COMSOL Multiphysics offers a range of solvers to solve various types of problems in physics.The choice of solver depends on the type of physics being modeled, the complexity of the problem, the sought-for accuracy, and the available computational resources.To solve the equations of the two-fluid turbulence model, a fully coupled approach was used with the direct solver algorithm (PARDISO).Newton's iterative method with a damping factor of 0.1 was used.The iterative process for the problem of flow around a flat plate with zero pressure gradient lasted up to 250 iterations, and for the remaining problems, the iterative process continued up to 350 iterations.The tolerance factor is 1, the residual factor is 1000.
The standard SST turbulence model uses standard COMSOL Multiphysics solvers.The following boundary conditions were set for the SST model: where L is the characteristic size of the streamlined body.The remaining boundary conditions were specified in the standard way.
For the problem of an aerodynamic profile, the initial distribution of velocity and pressure was given by the potential field of velocities.Assuming an irrotational inviscid flow, ϕ the velocity potential is defined as The velocity potential must satisfy the continuity equation for an incompressible flow ∇u = 0 .The continuity equation can be represented as the Laplace equation which is the potential flow equation.
After calculating the velocity potential, the pressure can be approximated using the Bernoulli equation for stationary flows:

Flat plate with zero pressure gradient
The main purpose of this experiment is to test the implementation of the two-fluid turbulence model in Comsol Multiphysics and compare the obtained results with the experimental data presented on the NASA website 7 .In the calculation, the plate length was L = 2 m with Reynolds number Re = 5,000,000 per unit length.In this case, the maximum thickness of the boundary layer is approximately 0.03 L, so the computational grid height was removed by distance y = L, which is sufficient to have little effect on the results.The turbulence intensity of the oncoming flow was 0.039%.The boundary conditions are shown in Fig. 1a, and an illustration of the computational grid and domain are shown in Fig. 1b.A computational grid of 69 × 49 in size, presented on the NASA website was used 7 .
Below are comparisons of the obtained numerical results with known experimental data.Figure 2 shows: a) the dependence of the friction coefficient along the plate; b) the dimensionless longitudinal flow velocity as a function of the dimensionless distance to the plate, as well as the results of Cole's theorie 40,41 .
Here C f is the coefficient of friction of the plate: The solid line in Fig. 3 shows the results of numerical calculation for the dimensionless longitudinal flow velocity depending on the dimensionless distance to the plate.Dimensionless velocities and distance were determined by the following formulas:  Figures 2 and 3 show that, for all Reynolds numbers, the model well describes both laminar and turbulent zones 22 .

Airfoil DSMA661
The DSMA661 airfoil is designed for low Reynolds flow.It was developed by the Delft University of Technology in the Netherlands.The DSMA 661 airfoil is commonly used in unmanned aerial vehicles (UAVs), especially those operating at low speeds.This airfoil has a relatively large maximum thickness of 16% of the chord length.It has moderate camber and is designed to provide good lift and low drag at low Reynolds numbers.The DSMA 661 airfoil is known for its stable and predictable behavior, making it suitable for a variety of UAV missions.
The DSMA661 airfoil 42 provides another opportunity to evaluate the implementation of the proposed twofluid model in Comsol Multiphysics.Four different two-dimensional grids were used in the work.Each coarser grid represents exactly every second point of the next finer grid, ranging from the 1121 × 193 fine grid to the coarsest 141 × 25 grid, presented on the NASA website 7 .There are 305 points on the fine grid along the trace from the trailing edge of the profile to the outflow boundary (39 points on the coarsest grid).The main results were taken for a grid of 561 × 97, and the remaining computational grids were used to check the grid convergence.The boundary conditions are shown in Fig. 4a, and an illustration of the computational grid and computational domain is shown in Fig. 4b.
For the problem at hand, the results of the well-known SST turbulence model are also obtained.The results are obtained with a free stream flow velocity of 18 m/s at zero angle of attack, which corresponds to a Reynolds number based on a chord of 1.2 million.For both models, the free flow turbulence intensity was 0.088%.
The distribution of the surface pressure coefficient on an airfoil is characterized by a change in pressure on its surface depending on the distance from a certain point.Typically, the analysis uses the surface pressure coefficient Cp, defined as the ratio of the pressure difference between a point on the surface of the profile and the pressure of the free flow to the dynamic pressure of the free flow.
where p is the pressure at a point on the profile surface, P ∞ is the pressure of the free flow, ρ is the density of the free flow, U 0 is the velocity of the free flow.
The distribution of the surface pressure coefficient on an airfoil can be used to analyze its aerodynamic characteristics such as lift, drag coefficient, etc. Figure 5a shows the distribution of the surface pressure coefficient C p of the DSMA661 profile at an angle of attack α = 0 0 .
The distribution of the coefficient of skin friction on an airfoil is characterized by a change in the friction force on its surface depending on the distance from a certain point.The coefficient of skin friction C f is defined as the ratio of the friction force acting on the surface of the profile to the dynamic pressure of the free flow.
where F is the friction force acting on the profile surface, S is the profile surface area oriented along the flow.Figure 5b illustrates the numerical and experimental results for C f for both the top and bottom profile surfaces.
Figures 5 show that the results of both turbulence models practically coincide and are in good agreement with the experimental data.( 16) Table 1 presents the errors in the deviation of numerical results from experimental data for C f and C p .Figures 6 and 7 show the longitudinal velocity U/U 0 profiles along the top (Fig. 6) and along the bottom (Fig. 7) surfaces of the profile at different sections along the flow.
As can be seen from Figs. 6 and 7, the results of both models are close to the experimental results.Below are the numerical results of the models and experiment for the profiles of longitudinal velocity U/U 0 (Fig. 8) and turbulent stresses u ′ ϑ ′ /U 2 0 (Fig. 9) in the wake after the profile at sections x/c = 1.01, 1.05, 1.20, 1.40, 1.80 and 2.19.
Table 2 presents the errors in the deviation of the calculated results from the experimental data in the cross section x/c = 2.19 for U/U 0 .
Figure 8 shows that the results of both models coincide with each other and with the experimental results for the average computational grid 561 × 97 with high accuracy.
For turbulent stress u ′ ϑ ′ /U 2 0 , there is also a good agreement between the results of models and experiment.Figure 10 shows the contours for longitudinal velocity and turbulent stresses.
To check the grid convergence, Fig. 11 shows the numerical results when the computational grid is changed.The results in section x/c = 0.99 for the bottom surface of the profile are shown.www.nature.com/scientificreports/In Fig. 11 shows that the results of both models do not depend on the resolution of the computational grid.
Figure 12 shows the change in the lift coefficient C L from the resolution of the computational grid.
It is also clear from Fig. 12 that the results of both models are less sensitive to the resolution of the computational grid.

Airfoil NACA 4412
The NACA 4412 airfoil is an airfoil with a maximum thickness of 12% of the chord length, located at 40% of the chord length.It is widely used in various fields, including aircraft wings and propeller blades.This section presents the validation of a two-fluid turbulence model for the NACA 4412 airfoil 43,44 .The case is considered for the Mach number M = 0.09, the angle of attack α = 13.87 • and the Reynolds number Re = 1,520,000.449 × 129 grid is presented in 7 .The boundary conditions are shown in Fig. 13a, and an illustration of the computational grid and domain is shown in Fig. 13b.
Figure 14 shows the isolines of the flow velocity.
It can be seen from this figure that both models show close results to the experimental data.
Figure 15 shows the distribution of the surface pressure coefficient C p on the NACA 4412 airfoil.From this figure it can be seen that the results of the two-fluid model are in better agreement with the experimental data.This can be explained by the fact that near the edge of the profile as shown in Fig. 14, a recirculating flow movement occurs, which causes anisotropic turbulence.The SST model uses the Boussinesq hypothesis, which is valid for isotropic turbulence.Therefore, under anisotropic turbulence, the SST results are somewhat worse.As for the two-fluid turbulence model, as shown in previous works, it is capable of describing anisotropic turbulence with great accuracy.
Table 3 presents the deviation of model results from experimental data for C p .The controls of the airliners are mainly located on the edge of the profile.Therefore, the accuracy of calculating the distribution of the pressure coefficient near the profile edge is of great practical importance.In this regard, Fig. 16 shows the distribution of the pressure coefficient at trailing of the edge.From this figure it can be seen that the results of the two-fluid model are significantly better than the results of the SST model.
Figure 17 shows the same results using different models for the NACA4412 airfoil at an angle of attack of 12 0 and Reynolds number Re = 1.64•10 6 .These results are obtained from the work 45 .
The results presented in Figs.16 and 17 show that the accuracy of the two-fluid model is higher than that of the RANS models and no worse than that of the LES and DES models.www.nature.com/scientificreports/ Figure 18 shows the longitudinal velocity U/U 0 profiles along the upper surface of the NACA 4412 airfoil at different sections downstream.
Here, too, there is a correspondence between the results of the models and experimental measurements.It can be seen from the figure that the velocity profiles according to the results of the SST model slightly deviate  4 presents the relative deviations of the model results from the experimental data in the cross section x/c = 0.6753 for U/U 0 .
Figure 19 shows the turbulent stress u ′ ϑ ′ /U 2 0 along the top surface for the NACA 4412 airfoil in different sections.
From Fig. 19 it is clear that the turbulent stresses correspond to the experimental data worse than the averaged values.This is due to the fact that they have small values and to improve compliance it is necessary to use calculation schemes of a higher order of accuracy.
Calculation time for problem NACA 4412: 420 iterations and 1740s in the SST model, and 380 iterations and 1588 s in the two-fluid model.All calculations were performed on a computer with a 2.5 GHz quad-core Intel i5-7300HQ processor, 16 GB DDR3 RAM, a 1024 GB hard drive, and Windows 7 (64-bit).

Conclusion
This article demonstrates the capability of the two-fluid turbulence model in the Comsol Multiphysics software package, which uses the finite element method.To stabilize the model equations, a stabilizer based on the Galerkin least squares method was used.For validation of the model and verification of the numerical algorithm, the problems of flow past a flat plate with a zero pressure gradient, as well as DSMA661 and NACA 4412 airfoils with angles of attack of 0 and 13.87 degrees, respectively, are considered.For the first time in the Comsol Multiphysics program, a two-fluid model was introduced and the results were obtained.To verify the created program, the results obtained are compared with the results of other models, as well as with experimental data.It is shown that the results of the two-fluid model correspond better to experimental data than the results of other RANS models and are no worse than the results of the LES and DES models.In addition, it is shown that the two-fluid model requires less computational resources than the SST model.Therefore, the two-fluid model can be recommended for calculating engineering problems of turbulent hydrodynamics.

Figure 3
Figure 3 shows the profiles of the dimensionless longitudinal velocity at different Reynolds numbers in two sections: (a) x = 0.97 m and (b) x = 1.97 m.The solid line in Fig.3shows the results of numerical calculation for the dimensionless longitudinal flow velocity depending on the dimensionless distance to the plate.Dimensionless velocities and distance were determined by the following formulas:

Figure 1 .
Figure 1.Flat plate with zero pressure gradient (a) boundary conditions and (b) computational grid and domain.

Figure 2 .
Figure 2. Dependence of the coefficient of friction along the plate (a), profile of the dimensionless longitudinal flow velocity on the dimensionless distance to the plate (b).

Figure 11 .
Figure 11.Checking the grid convergence of turbulence models.

Figure 12 .
Figure 12.Change in the lift force coefficient C L from the resolution of the computational grid.

Figure 15 .
Figure 15.Results for surface pressure coefficient.

Figure 16 .
Figure 16.Pressure coefficient distribution at trailing of the edge.

Figure 17 .
Figure 17.Pressure coefficient distribution near the profile edge.

Table 1 .
Advances in LES of Complex Flows.Here δ is the relative deviation.

Table 2 .
Deviation of tu.Here δ is the relative deviation.

Table 3 .
The relative.Hear δ = C p.exp − C p.model .C p upper wall

Table 4 .
The Cp model.Here δ is the relative deviation.